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The  paper  addresses  the  question  of  the  reliability  of  engineering  compu¬ 
tations.  It  brings  a  set  of  paradoxical,  unexpected  results  which  shows  that 
the  common  practice  can  lead  to  unreliable  results  and  conclusions.  The 
theory  and  implementation  of  the  analysis  of  elasticity  problems  with 
stochastic  input  data  (loads,  domain,  coefficients)  are  outlined.  Numerical 
examples  Illustrate  the  ideas  and  results. 
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ABSTRACT 

The  paper  addresses  the  question  of  the  reliability  of  engineering 
computations.  It  brings  a  set  of  paradoxical,  unexpected  results  which  shows 
that  the  conmon  practice  can  lead  to  unreliable  results  and  conclusions. 

The  theory  and  implementation  of  the  analysis  of  elasticity  problems  with 
stochastic  input  data  (loads,  domain,  coefficients)  are  outlined.  Numerical 
examples  illustrate  the  ideas  and  results. 
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I .  INTRODUCTION 

/  Shape  optimlzaclon  In  the  'structure  mechanics  became  to  be  in  the  center 


of  the  research  and  applications.  Many  papers  and  books  dealing  with  this 
subject  appeared  and  special  conferences  dealt  with  these  problems.  The 
research  is  directed  toward  theoretical  questions  as  existence  and  characteri¬ 


zation  of  the  optimal  design,  bounds  for  the  optimized  values,  numerical 
treatment  of  the  optimal  design  problems,  etc.  We  will  not  mention  here  the 

y 

recent  vast  available  literature.  We  mantldn  only  [10]  [13]  [18]  [21]  as 


m 


examples. 


We  will  address  in  this  papery|the  problem  of  th^reliablllty  of  the 
conclusions  based  on  the 'computational  analysis  and  their  relation  to  the 


problems  of  the  optimal  design. 

By  reliabyjjty-we  mean  here  that  the  conclusions  sufficiently  accurately 
describe  the  physical  reality. 


The  problem  of  optimal  design  consists— in  principle — in  the  comparison 
of  the  solutions  of  states  from  the  set  S  of  admissible  states  (for  example. 


solution  of  the  problems  from  the  set  of  admissible  domains)  and  the  selection 


of  Che  "optimal"  state  (e.g.  the  domain)  for  the  engineering  design.  It  is 


obvious  that  such  a  selection  can  be  successful  only  if  the  solution  of  ever^ 


state  is  uniformly  reliable  with  respect  to  the  entire  set  S.  This  require¬ 


•  .  •  J 

m 
.  • .  “ .  <'1 
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ment  creates  a  serious  difficulty  because  we  are  used  to  solve  in  practice 
numerically  the  simplified  mathematical  formulation  of  the  problem  and  have 


experience  only  with  small  limited  set  of  practical  problems.  Hence,  It  is 
essential  to  analyze  and  to  be  explicitly  aware  of  the  assumptions  used  in  the 
derivation  of  the  model  and  its  numerical  treatment  and  the  limitations  we 
have  to  deal  with. 


i 
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From  what  we  said ,  it  is  obvious ,  that  we  have  to  focus  on  the 
reliability  of  the  analysis  of  the  single  states  and  its  uniformity  over  the 
(entire)  set  S  of  admissible  states.  This  has  always  to  be  the  starting 
point  of  the  assesment  of  the  validity  of  the  optimal  design. 

The  reliability  of  the  computational  treatment  of  the  single  states 
depends  on 

a)  mathematical  model 

b)  reliability  of  the  input  data 

c)  reliability  of  the  numerical  treatment. 

These  three  aspects  are  of  course  closely  related. 

In  this  paper  we  will  address  the  questions  not  in  general ,  but  on  few 
concrete  engineering  examples.  Ito  restrict  ourselves  to  the  examples  which 
are  relatively  simple  from  the  engineering  (although  not  mathematical)  point 
of  view,  to  present  the  ideas  in  a  most  clear  way. 

2.  THE  PROBLEM  OF  A  TUBE  WITH  A  STIFFENED  SURFACE 

Let  us  consider  the  problem  of  a  tube  (plane  strain)  with  stiffened 
outer  surface.  We  shall  assume  that  the  stiffener  is  bending  free.  We  are  in 
general  interested  in  the  design  resp.  the  influence  of  the  changes  of  the 
outer  surface  Tg.  The  scheme  of  the  problem  is  shown  in  Fig.  2.1. 

We  will  formulate  the  problem  as  linear  elasticity  ( plane  strain) 
problem.  The  formulation  is  the  standard  one,  based  on  the  minimization  of 
Che  (cummulative)  energy  of  Che  tube  and  the  (tension)  energy  of  the  rein¬ 
forcement  on  Tq. 

2.1.  The  reliability  of  the  mathematical  model 

We  will  analyze  the  case  when  the  outer  boundary  Tg  is  a  regular  m- 
polygon  and  T,  is  a  concentric  circle.  See  Fig.  2.2.  The  domain  is  denoted 


by  The  circular  domain  JIq  (see  Fig.  2.3)  is  obviously  the  limiting 

case  when  m  ->■  “>.  Assuming  that  the  unit  hydrostatic  pressure  is  given  on  the 
inner  boundary  and  that  the  stiffener  is  infinitely  rigid  in  tension,  the 
problem  of  linear  elasticity  can  be  formulated  on  one  sector  only  with  the 
boundary  conditions  shown  in  Fig.  2.4.  On  the  sides  A-B  and  C-D  the  symmetry 
conditions  are  prescribed.  On  the  side  A-D  the  boundary  conditions  describe 
the  behaviour  of  the  stiffened  side  and  on  the  side  B-C  the  tractions  are 
prescribed . 

Denote  by  resp.  (uq.Vq)  the  solution  (displacement)  on 

resp.  Qq. 

From  Che  physical  grounds  we  have  Co  expect  that 
as  m  +  If  (u^.Vjjj)  ^  (uq,Vq),  we  have  to  have  doubts  about  the  relia¬ 
bility  of  the  model. 

We  have 

THEOREM  2.1. 

(2.1)  lira  (u  ,u  )  =»  (u  ,v  )  /  (up^,v_). 

m  m  oo'  00  U  u 

m-«® 


For  the  analysis  of  the  dependence  of  the  solution  on  small  changes  of  the 
domains  we  refer  e.g.  to  [1]  [3]  [4], 

Theorem  2.1  shows  that  Che  used  model  of  linear  elasticity  is  unreliable 
at  least  if  the  outer  boundary  is  not  smooth  and  hence  it  practically  cannot 
be  used  for  optimal  design  when  the  admissible  domains  have  not  sufficiently 
smooth  boundary.  We  see  here  that  the  uniform  reliability  (with  respect  to 
m)  is  clearly  violated. 


The  limiting  solutions  Cu^,v^)  and  found 


THEOREM  2.2 

.  The  solutions  (u  ,v  )  and 

(uq.vq)  are  radially  symmetric 

Denoting  by 

a  resp.  a„  the  stresses 

r  0 

in  polar  coordinates  we  get: 

(2.2) 

A  .  -n 

Or  =  “2 

"e  ■  -  \  » 

r 

r 

with 

(2.3) 

B  *  0 

aa 

(2.4) 

^  _  (l-v)a^b^ 

(l+v)a^ 

^  (l-v)b^  +  (l+v)a^’ 

2  2  * 
(1-v)d  +  (l+v)a^ 

Table  2.1  gives  the 
for  the  solution  (u  ,v  ) 

40  00 

0.3.  We  see  clearly  that 
tially  different. 


values  of  the  stresses  a^.  and  Oq  on  the  line  C-D 
and  (uq,Vq)  when  a=0.3,  b=1.0  and  v  » 
the  solutions  (u^,v^)  and  ('^»Vq)  are  essen- 


TABLE  2.1 


r 

(u«,v, 

-) 

c 

o 

)> 

‘^r 

<^0 

0.3 

-1.000 

1.000 

-1.000 

0.7135 

0.4 

-0.5635 

0.5625 

-0.6251 

0.3387 

0.5 

-0.3600 

0.3600 

-0.4516 

0.1652 

0.6 

-0.2500 

0.2500 

-0.3574 

0.0710 

0.7 

-0.1836 

0.1836 

-0.3005 

-0.0142 

0.8 

-0.1406 

0.1406 

-0.2637 

-0.0227 

0.9 

-0.1111 

0.1111 

-0.2384 

-0.0480 

1.0 

-0.0900 

0.0900 

-0.2208 

-0.0609 

Before  discussing  the  probable  reason  for  this  paradox,  let  us  mention 


THEOREM  2.3,  Let  be  the  solution  on  0^^  (i.e.,  m-sided  polygon), 

ra  >  A.  Then  the  solution  has  a  singularity  in  the  neighborhood  of  the  point  A 
(vertex  of  the  polygon)  (see  Fig.  2.2  and  2.4)  and 


(2.5) 


u 

m 

V 

m 


+  higher  order  terms 


where  (r,0)  are  Che  polar  coordinates  with  the  origin  In  A  and 

are  smooth  functions  in  0.  H 

Theorem  2.3  shows  that  the  solution  has  a  strong  singularity  in  the 
neighborhood  of  A  and  the  strains  and  stresses  are  there  unbounded.  This 
obviously  violates  the  basic  assumptions  of  the  linear  elasticity  model  and 
has  unexpected  consequences. 

We  are  making  the  following  conjecture  (unproven): 

If  (u  ,v  )  is  the  solution  of  a  nonlinear  problem,  then  lim(u  ,v  )  •► 
m  m  . .  ra  m 

^  ’^0  ^  * 

This  leads  to  the  following  conclusion.  If  the  set  of  admissible 
domains  has  unsmooth  outer  boundary,  Chen  it  is  necessary  to  use  nonlinear 
theory  of  elasticity  in  the  optimal  design  problems.  The  linear  elasticity 
leads  to  unreliable  results  and  conclusions. 


2.2.  The  reliability  of  the  numerical  solution 

As  we  have  seen  in  Theorem  2.3,  the  solution  has  a  very  strong  singu¬ 
larity  (note  that  in  the  case  of  a  crack  the  singularity  of  Che  solution  is 
r'^2  )  which  makes  computation  very  difficult  for  larger  m.  Tlie  computation 
we  present  has  been  make  by  the  code  PROBE  which  uses  v  and  h-p  version  of 
Che  finite  element  method.  See  [23]  (24],  For  the  theoretical  aspects  we 
refer  to  [6]  [7]  [14],  The  mesh  has  to  be  strongly  refined  in  the  area  of  the 


singularity  if  reliable  results  have  to  be  obtained.  See  [14]  [22].  For 
the  p-version  (i.e.  when  there  is  no  strong  refinement  at  A)  the  energy  norm 
of  the  error  is  Bell  »  Cp  See  [6]  [9].  For  the  h-version  without 

properly  refined  mesh,  the  situation  is  still  worse.  For  the  properly  refined 

mesh  the  rate  of  convergence  in  the  first  phase  (p  not  large)  is  exponential. 
See  [14],  The  mesh  we  used  is  shown  in  Fig.  2.5  (a  =0.3,  b  =  l,  v  =  0.3). 

In  Table  2.2  we  show  the  stresses  in  the  points  P,  Q,  R,  S  (see  Fig. 

2.5)  for  ra  =  8,  16,  32  for  various  degrees  p  of  elements.  We  see  that  the 

solution  is  close  to  the  limiting  value  of  m  =  ».  Although  a  =  = 

0  we  see  that  a  deteriorates  from  m  =  16  to  m  =  32.  The  reason  is  that 
the  quality  of  the  numerical  solution  deteriorates  with  m  -►  ®.  This  deteri¬ 
oration  is,  for  example,  visible  from  the  Table  2.3  where  the  computed  strain 
energy  for  various  p  and  m  is  given.  We  see  clearly  a  much  larger  change 
in  the  energy  for  m  =  32  than  for  m  =  8  when  increasing  the  degree  p. 

This  indicates  much  larger  error  for  m  =  32  than  for  m  =  8.  The  strength 
of  the  singularity  is  j-2/(m-2)  which  is  so  strong  that  without  special  care 

no  reasonable  accuracy  can  be  achieved  for  ra  =  32.  For  ra  =  16  the  error  in 

the  energy  norm  is  expected  to  be  2-4%,  and  7-10%  for  ra  =  32  in  our 
computations . 

Table  2.4  shows  the  values  of  the  maximal  principle  stress  in  the 
points  B-F  and  B-F  (see  Fig.  2.5).  We  clearly  see  that  the  stresses  are 
very  large  in  the  neighborhood  of  the  vertices  and  with  m  +  “  the  stresses 
are  increasing  (because  the  strength  of  the  singularity  is  increasing).  This 
also  clearly  indicates  the  likely  reason  for  the  paradox  we  mentioned  above. 

We  see  that  not  only  the  mathematical  model  but  also  the  quality  numer¬ 
ical  solution  is  very  nonuniform  with  respect  to  small  changes  of  boundary 
(which  does  not  have  sufficient  smoothness). 


o  ^ 


TABLE  2.4 


1 

P 

Maximal  principal  stress 

m  =  8 

m  =  16 

m  =  32 

F 

8 

-0.3807+0 

-0.3758+0 

-0.3624+0 

7 

-0.3814+0 

-0.3767+0 

-0.3670+0 

6 

-0.3802+0 

-0.3716+0 

-0.3494+0 

E 

8 

-0.1360+1 

-0.1907+1 

-0.2043+1 

7 

-0.1358+1 

-0.1916+1 

-0.2076+1 

6 

-0.1355+1 

-0.1879+1 

-0.1990+1 

D 

8 

-0.4805+1 

-0.9498+1 

-0.1150+2 

7 

-0.4810+1 

-0.9502+1 

-0.1153+2 

6 

-0.4748+1 

-0.9394+1 

-0.1119+2 

C 

8 

-0.1709+2 

-0.4695+2 

-0.6393+2 

7 

-0.1709+2 

-0.4683+2 

-0.6383+2 

6 

-0.1769+2 

-0.4585+2 

-0.6176+2 

B 

8 

-0.5296+2 

-0.1764+3 

-0.2722+3 

7 

-0.6709+2 

-0.2762+3 

-0.4340+3 

6 

-0.4972+2 

-0.1442+3 

-0.1931+3 

F 

8 

-0.3851+0 

-0.4044+0 

-0.5062+2 

7 

-0.3855+0 

-0.4048+0 

-0.4998+0 

6 

-0,3821+0 

-0.3947+0 

-0.4767+0 

E 

8 

-0.1356+1 

-0.2111+1 

-0.3358+1 

7 

-0.1358+1 

-0.2135+1 

-0.3427+1 

6 

-0.1335+1 

-0.2083+1 

-0.3351+1 

D 

8 

-0.4818+1 

-0.1130+2 

-0.2136+2 

7 

-0.4825+1 

-0.1147+2 

n  o  1  o  / 

“vj  •  OH-r^ 

A 

-0.4748+1 

-0.1121+2 

-0.2136+2 

C 

8 

-0.1740+2 

-0.6362+2 

-0.1400+3 

7 

-0.1747+2 

-0.6459+2 

-0.1435+3 

6 

-0.1722+2 

-0.6349+2 

-0.1408+3 

B 

8 

-0.6361+2 

-0.3590+3 

-0.9021+3 

7 

-0.6481+2 

-0.3824+3 

-0.9606+3 

6 

-0.6437+2 

-0.3636+3 

-0.8924+3 

3.  THE  PROBLEM  OF  THE  PLATES  AND  SHELLS 

In  Section  2  we  addressed  the  problem  of  the  reliability  of  the  linear 
elasticity  model.  Models  of  plates  and  shells  are  two  dimensional  although 
obviously  the  original  problem  is  three  dimensional.  Hence,  we  will  assume 
that  the  three-dimensional  linear  elasticity  formulation  is  reliable  and  will 
analyze  only  the  effects  of  the  dimensional  reduction  from  3  to  2  dimensions 
and  the  implication  for  the  optimal  design. 

3.1.  The  problem  of  the  simply  supported  plate 

Let  us,  for  simplicity,  assume  that  we  are  concerned  with  the  case  when 
the  Poisson  ratio  v  =  Q.  The  plate  problem  (with  uniform  thickness  h)  can  be 
formulated  in  various  ways.  Let  us  mention  the  "projection  method"  when  we 
assume  an  "ansatz"  and  use  it  in  the  variational  principle  by  minimizing  the 
energy.  This  approach  is  sometimes  called  Kantorowich  method  (see  [15]). 
Denoting  u,  v,  w  the  displacement  components,  we  shall  consider  two 
"ansatzes" ; 

1)  The  K  (Kirchhoff)  model 
(3.1a)  u(x,y,z)  =  ~z  (x,y)  , 

(3.1b)  v(x,y,z)  =  -z  (x,y)  , 

(3.1c)  w(x,y,z)  =  w(x,y). 

Using  this  ansatz  in  the  potential  energy  principle  we  get  the  usual 
formulation 


(3.2) 


El  AA  w  =  f 


9 


and  the  simple  support  is  obtained  by  minimization  of  the  energy  with  the  onl^ 
constraint  w  =  0  on  the  boundary  T  of  the  domain. 

2)  The  R-M  (Reissner-Mindlin)  model 


(3.3a) 

u(x,y,z) 

=  -z<p(x,y) 

(3.3b) 

v(x,y,z) 

=  -zt|)(x,y) 

(3.3c) 

w(x,y,z) 

=  w(x,y) 

V\-\* 


Utilizing  (3.3)  in  the  expression  for  the  three-dimensional  potential  energy 
and  imposing  the  (only)  constraint  w  =  0  at  T  we  obtain  a  system  of  three 
differential  equations  of  second  order  in  contrast  to  one  equation  of  fourth 
order  in  the  K-raodel. 

The  dimensional  reduction  has  been  analyzed  in  the  asymptotic  way  when 
h  >  0  and  the  solution  is  smooth.  See  e.g.  [11]  [12]  [20].  In  this 
asymptotic  frame  we  cannot  distinguish  between  the  two  mentioned  models. 
Physically  the  R-M  model  is  taking  into  account  the  shear  stresses  while  the  K 
model  neglects  them.  Let  us  once  more  assume  that  is  the  regular  ra- 

polygon  inscribed  in  the  circle  of  radius  a,  Uq  be  the  circle  with  radius 
a  and  let  us  consider  the  problem  of  uniformly  loaded  (by  load  p)  simply 
supported  plate. 

Denote  by  w^^,  resp.  and  t resp.  ^  the  solution 


^  m  u  m  m  m 

of  the  K  and  K"M  model  on  and  Qq.  Then  we  have 


THEOREM  3.1. 
(3.4a) 

(3.4b) 


”m  ^  '^0 


10 


I  We  can  compute  the  limiting  solution  w„  and  ( analytically.  For 

the  K-model  we  have 

I  wjO.O)  -  - 

(3.5b)  Wq(0,0)  *  ~  ZU  ^  ^ 

I 

and  for  the  R-M  model  we  have 

(3.5c)  w„(0.0)  =  Wq(0,0)  = 

! 

where  1  =  h'’/12,  F  =  h  are  the  moment  of  inertia  and  the  thickness, 
respectively,  and  E  is  the  modul  of  elasticity. 

I  Theorem  3.1  and  the  formulae  3.5a-3.5c  show  that  the  effects  of  the 

shear  stress  in  the  neighborhood  of  the  corners  are  essential.  Although  we 
discussed  only  the  problem  of  the  polygon  plate,  the  analysis  (see  [3]  (5]) 
covers  much  more  general  solutions  and  clearly  point  to  the  following 
conclussion: 

The  optimal  design  of  a  plate  has  to  be  based  on  the  R-M  and  not  the 
K-model . 

We  will  not  discuss  here  the  reliability  of  the  numerical  treatment.  Similar 
but  more  complicated  situation  occurs  in  the  case  of  the  shells. 

3.2.  The  problem  of  the  plate  with  a  variable  thickness 

Let  us  consider  a  plate  with  variable  thickness.  If  the  thickness  is 


very  slowly  varying  with  respect  to  the  thickness  of  the  plate,  then  the 
derivation  (dimensional  reduction)  can  be  made  in  the  same  way  as  for  the 


constant  thickness.  Nevertheless,  if  the  thickness  is  varying  rapidly,  then 
the  classical  derivation  is  not  valid.  Recently  a  theory  has  been  developed 
(see  [16]  [17])  which  shows  Important  relation  between  the  thickness  and 
thickness  variation  and  which  strongly  influences  the  reliability  of  the 
mathematical  model.  We  will  show  it  in  the  most  simple  setting.  Consider  the 
stiffened  plate  shown  in  Fig.  3.1.  The  main  idea  of  the  classical  plate 
derivation  is  to  consider  the  limiting  process  e  -►  0  and  apply  the  results 
for  e  >  0. 

X  <  X2 

We  can  assume  that  a  =  C^e  ,  b  »  consider  the  limiting 

process  e  ->■  0.  In  [17],  X^  =  X^  =  X  is  assumed  and  it  is  shown  that  we  get 
different  model  for  X  <  1,  X  =  1,  and  X  >  1. 

In  the  case  of  X  <  1  the  stiff ners  are  far  apart  when  e  -*■  0,  in  the 
case  X  >  1  they  are  close  together.  In  all  three  cases  the  dimensional 
reduction  leads  to  the  plate  formulation  with  effective  coefficients  depending 
on  the  value  of  X.  This  example  shows  that  optimal  design  based  on  one 
model,  say,  X  <  1  for  fixed  but  small  thickness  can  lead  to  a  design  when 
the  model  is  not  valid  (reliable)  anymore.  Using  a  proper  model  for  this 
design  and  redoing  the  optimal  design  once  more,  we  can  once  more  get  out  of 
the  range  of  the  reliability  of  the  model.  Hence,  we  have  to  consider  here 
simultaneous  design  optimization  and  the  model  selection.  For  important 
aspects  of  this  problem  directly  related  to  the  optimal  design,  we  refer  to 
[17]. 

4.  THE  PROBLEM  OF  A  SUPPORTED  CONSTRUCTION 

Let  us  consider  the  optimal  design  of  a  supported  construction  (see  Fig. 


4.1).  The  problem  is  how  to  model  the  support  in  the  point  B.  To  show  the 
difficulty,  let  us  consider  the  problem  shown  in  Fig.  4.2  and  solve  the  linear 


elasticity  (plane  stress,  v  »  0.3)  problem.  The  standard  finite  element 
modeling  is  to  make  constraint  v  =  0  at  the  node  located  in  the  support. 

This  modeling  is  incorrect  because  the  solution  strongly  depends  on  the  finite 
element  method. 

Let  M  be  the  moment  at  the  side  A-A'  and  Mjj  is  the  moment  computed 
by  the  finite  element  method.  Assume  that  the  size  h  of  the  maximal  element 
hj^  =  max  h  -►  0  as  N  ■►  *.  We  have 

THEOREM  4.1. 

(4.1)  lim  Mjj  =  Mq 

N-H» 

where  Mq  is  the  moment  when  there  is  no  support  (and  hence  Mq  can  be 
analytically  computed) . 

Theorem  4.1  shows  Chat  by  selecting  different  meshes  we  can  get 
completely  different  results  and  hence  optimal  design  will  strongly  depend  on 
the  used  mesh.  In  fact,  the  situation  is  still  more  complicated  because  M^ 

Mq  slowly  and  we  have  no  means  to  establish  how  reliable  the  solution  is. 

Before  discussing  this  effect,  let  us  show  the  computation  by  the  code 
PROBE.  The  used  mesh  is  shown  in  Fig.  4.3.  There  is  refinement  in  the  neigh¬ 
borhood  of  A-A'  and  especially  strong  refinements  is  in  the  neighborhood  of  B 
(see  Fig.  4.4).  We  did  use  two  meshes  A^,  with  smallest  ring  of  Che  radius 
a^,  and  A^  with  the  radius  a^.  Table  4.1  shows  the  moments  on  A-A'  and 
Table  4.2  shows  the  displacement  v  in  C.  Although  the  moments  and  the  dis¬ 
placement  are  significantly  smaller  than  that  of  the  unsupported  beam,  the 
mesh  dependence  is  obvious.  Note  that  Che  difference  between  the  values  ob¬ 
tained  by  the  mesh  A^  and  A^  is  nearly  independent  of  p.  The  reason  for 
the  effects  we  have  shown  is  that  the  support  is  not  correctly  modeled.  The 


reaction  is  a  point  force  which  leads  to  the  infinite  energy  and  a  infinite 
displacement  in  the  point  of  the  reaction.  The  infinite  displacement  at  the 
reaction  point  can  be  seen  from  the  analytical  solution  on  half  plane  with 
concentrated  load.  Hence,  reaction  has  to  be  zero  and  we  obtain  the  solution 
of  an  unsupported  beam. 


TABLE  4.1 


In  Table  4.3  we  show  the  displacements  in  the  points  and  B^^ 

computed  by  the  raesh  and  A5.  Realizing  that  the  distance  between  B 

and  B5  is  0.75  10”^  and  the  constraint  in  B  is  v  =  0,  we  see  niiraeri- 
cally  the  effect  mentioned  above.  This  clearly  shows  that  the  mathematical 
model  of  the  supported  beam  is  unrealiable  because  it  does  not  distinguish 
between  supported  and  unsupported  beam.  Hence,  a  more  sophisticated  model  of 


TABLE  4.3 


1 

P 

Mesh 

A4 

Mesh 

A5 

P 

0 

I 

N 

T 

Mesh 

A4 

Mesh 

^5 

®5 

8 

-  9.189 

®5 

-  9.188 

7 

-  8.777 

D 

-  8.777 

6 

-  8.313 

-  8.312 

5 

-  7.764 

-  7.763 

^4 

8 

-  9.231 

^4 

-  9.226 

7 

-  8.818 

-  9.831 

6 

-  8.352 

-11.49 

-  7.795 

-11.48 

5 

-  7.801 

-10.94 

-  7.123 

-10.94 

8 

msE/m 

h 

7 

6 

-11.57 

-14.64 

5 

-10.97 

-14.10 

8 

-15.72 

-15.48 

-18.61 

7 

-15.31 

-15.07 

-18.20 

6 

-14.85 

-14.05 

-17.73 

5 

-14.31 

-17.40 

-13.37 

-17.19 

8 

-19.60 

-22.52 

®i 

WS3M 

-21.27 

7 

-19.22 

-22.14 

-20.84 

6 

-18.78 

-21.17 

-20.36 

5 

-18.27 

-21.20 

-15.80 

-19.79 

the  support  is  needed.  Nevertheless,  we  will  not  discuss  here  the  question  of 
a  reliable  model.  (Usually  it  is  claimed  that  a  concrete  not  strongly  refined 
mesh  models  the  support.  It  is  obvious  that  without  a  reference  to  the  proper 
formulation  of  the  support  the  claim  has  no  firm  meaning.] 


5.  THE  PROBLEM  OF  THE  STOCHASTIC  INPUT  DATA 

The  basic  Input  data  describing  the  elasticity  problems  are:  the  do¬ 
main,  the  material  propertis  and  the  loads.  Assume  now  that  the  data  are 
stochastic  functions.  For  example,  the  boundary  of  the  domain  can  be 


■  vVc--.  L-.  -V 
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described  by  a  stochastic  function  which  expresses  the  uncertainty  of  fabri¬ 
cation.  Then  the  solution  is  also  a  stochastic  function.  In  addition,  the 
failure  crlterium  which  can  be  basis  for  the  optimal  design  is  always  a 
stochastic  one.  Hence,  we  have  combine  both  stochastic  characters  to  get 
desired  information.  Because  the  uncertainty  of  the  input  data,  the  disper¬ 
sion  of  the  results  can  be  significant.  Recently  we  developed  a  theory  of  the 
solution  with  stochastic  input  data  (see  [2]  (19]  and  forthcoming  papers)  and 
their  numerical  tratment  by  the  finite  element  method.  The  implementation  is 
based  on  the  code  PROBE  mentioned  earlier. 

5.1.  The  case  of  the  stochastic  load 

Let  us  consider  the  container  of  the  form  shown  in  Fig.  5.1.  The  side 


AB  is  loaded  by 

a  horizontal 

stochastic  function  = 

A(y  ,u>)  , 

0  <  y  <  H 

and 

we  will  assume 

that  K  (y.m) 

»  1  where  we  denoted  by 

X„(y,w) 

the  mean. 

The 

correlation  function  is  assumed  to  be 

(5.1)  K(ypy2)  =■  0.1^  e"®!  ^I'^a  I  * 

A  simulated  sample  of  the  load  from  a  given  probability  field  (o= 0.03781) 
is  shown  in  Fig.  5.2.  We  will  assume  that  the  linear  elasticity  provides 
reliable  results  for  all  loads  under  consideration.  The  solution  of  our  model 
problem  is  a  stochastic  function  with  the  mean  being  the  (deterministic) 
solution  for  the  mean  load. 

Our  aim  is  to  determine  the  variance  and  covariance  of  the  values  of 
Interest . 

Concerning  the  failure  criterion  we  will  assume  as  example: 

a)  The  criterion  of  stress  Intensity  factor  F  in  the  point  C  (see 


i  **,»■' 


‘V\ 


-  > 


Fig.  5.1) 


b)  The  failure  criterion  based  on  the  envelope  of  the  Mohr's  circles  in 
the  point  D  (see  Fig.  5.1). 

Knowing  the  stress  intensity  factor  as  random  variable  characterized  by 
its  mean  and  standard  deviation,  we  can  establish  the  probability  level  of  the 
failure  when  the  material  probabilistic  characterization  of  admissible  stress 
intensity  factor  is  given. 

The  criterion  b)  is  more  complicated.  We  need  here  the  correlation  of 
the  components  of  the  stress  tensor  which  allows  us  to  compute  not  only  the 
mean  Mohr's  circles  but  also  its  perturbation  in  every  point  of  the  circle 
which  for  a  given  probability  level  has  an  elliptic  character.  The  envelope 
of  these  ellipses  will  be  compared  with  the  admissible  failure  curve  (see  Fig. 
5.3) . 

The  concrete  computation  of  our  model  problem  has  been  computed  using 
the  program  PROBE  and  the  mesh  shown  in  Pig.  5.4.  (The  refinement  in  the 
neighborhood  of  the  reentrant  corners  is  not  shown.) 

In  Table  5.1  we  show  the  mean  values  and  the  standard  deviation  Sd(F) 
of  the  stress  intensity  factor  F  in  the  point  C  in  dependence  on  p. 

For  the  technique  used  in  PROBE  for  the  computation  of  the  stress 
intensity  factor,  see  [8]. 


TABLE  5.1 


-46.5958 

-51.7433 

-49.3796 

-49.0721 


sd(F) 


3.71714 

3.92931 

3.94039 

3.91575 


.••V- 
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The  Mohr'  circles  for  the  probability  level  90%  are  shown  in  Fig.  5.5 


5.2.  The  problem  of  stochastic  boundary 

The  problem  of  stochastic  boundary  is  more  complicated  but  it  can  be 
transformed  to  the  case  of  stochastic  load.  Let  us  consider — for  simplicity 
of  the  exposition — the  problem  of  a  symmetric,  cracked  panel  (plane  strain,  v 
=  0.3)  shown  in  Fig.  5.6  and  assume  that  the  deterministic  traction  T  at  the 
boundary  is  such  that  the  exact  stress  tensor  is  given  by  the  following 
formulae 


(5.2) 


=  (2rrr)  cos  y  (1  ~  sin  sin 

=  (2iTr)  cos  Y  (1  +  sin  y  sin  -j^) 

=  (2Trr)  sin  |-  cos  j  cos  —  . 


These  functions  are  symmetric  mode  functions  of  the  stress  intensity  factor. 
Let  us  assume  that  the  side  A  is  perturbed  by  a  stochastic  perturbation  so 
that  the  boundary  is  given  by  the  function  y  =  1  +  A(x,u)),  -1  <  x  <  1, 

A(±l,u)  =  0,  where  A(x,(i))  is  the  stochastic  function  with  the  correlation 
function  K2(xj,X2). 

We  use  in  our  model  problem 

|x,-x  1 

(5.3)  K]^(x^,X2)  =  f( - ^ - 

where 

f(?)  =■ 


X 

-  0.  <  1,  i  =  1,2, 

8-l5C^(2~0^ 

720 


A  simulated  sample  of  the  perturbance  is  shown  in  Fig.  5.7.  Our  aim  is  to 
find  the  stress  intensity  factor  F  and  its  standard  deviation  Sd( F)  caused 
by  the  random  boundary  and  the  stresses  and  their  variances  and  covariances 
in  (0.1,  0.9) . 

Before  addressing  this  problem,  we  have  to  know  how  the  traction  will 
change  when  the  domain  is  changing  so  that  the  equilibrium  is  always  guaran¬ 
teed.  We  will  assume  to  this  end  that  functions  aj^(x,y)  ,  T^y(x,y)  and 
?fy(x,y)  are  defined  in  the  neighborhood  of  the  side  A-B  such  that 


(5.4) 


0 


3o 

3y 


0 


and  a„  =•  and  on  A-B  where  (x  ,a  )  is  the  given  traction 

y  y  xy  xy  xy  y 

vector  at  A-B.  If  now  point  D  lies  on  the  perturbed  boundary,  then  the 
traction  vector  T  is 


(5.5) 


where  (np  02)  is  the  outer  normal  to  the  perturbed  boundary.  This 
guarantees  the  equilibrium  condition  for  every  perturbation. 

Assume  that  the  magnitude  of  the  perturbance  is  X.  Then  we  have 


THEOREM  5.1.  The  solution  of  (up  to  higher  order  terms  in  X)  the  perturbed 
problem  is  the  solution  of  the  original  domain  with  the  modified  load  T 


T»TTr 


(5.6)  T  =  Tq  -  A(x.ai)  ^ 


-T  +  T 
xy  xy 


-a  +  a 

y  y 


+  A''(x,u) 


a  +  a 

X  X 


-T  +  T 

xy  xy 


where  Tn  =  is  the  traction  on  A-B. 

U  xy’  y 


*  -•  V  *1 
s  V  V 


Theorem  5.1  gives  us  Immediately  the  possibility  to  solve  the  problem  in  the 
same  vein  as  in  the  previous  section. 

We  have  used  in  our  model  problem 


yx.y)  =  <T^(0,1) 


T  (x,y)  =  T  (x,l) 

xy  *•’  xy  * 


^y(x,y)  =  ffyCx,!)  -  (y-1)  (x,l). 


We  used  the  program  PROBE  and  for  p  =  8  we  obtained 

a)  The  stress  intensity  factor  F: 

the  mean  value  F  =  0.99830  (exact  value  F  =  1) 
the  standard  deviation  sd(F)  =  2.54(-4). 

b)  The  stress  in  the  point  A  =  (0.1,  0.9): 

the  mean  value:  a  =  0.1426,  a  =  0.4821,  t  =  0.1206 

X  y  ’  xy 

the  standard  deviation  sd(a„)  =  0.48418(-2),  sd(a  )=  0.35178(-2), 

X  y 


sd(T^y)  =  0.2088(-2) 

the  covariance 


c(CTx>‘^y>  =  0.1697 (-4),  c(a^,  x^^) 


0.9493(-5)  c(a  ,T  )  =  0.7019(-5) 
y  xy 

the  normalized  covariance  p(a  ,a  )  =  0.9966,  p(a  ,t  )  =  0.9389, 

x’  y  X  xy 

pCOy.x^y)  =  0.9555. 
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We  see  that  the  variance  of  F  is  much  smaller  than  the  variance  of  the 
stress  in  the  point  A.  If  the  failure  criterium  is  based  on  the  stress 
intensity  factor  F,  then  it  is  (in  our  case)  practicaly  uninfluenced  by  the 
uncertainty  of  the  boundary.  If  the  failure  creterium  is  based  on  the  Mohr 
circle  in  A,  then  it  is  much  more  sensitive  to  the  uncertainty  of  the 
boundary.  This  shows  very  clearly  that  the  same  uncertainties  can  lead  to  the 
uncertainties  of  different  magnitude  in  the  failure  criterium  parameters. 

Let  us  mention  that  we  need  in  (5.6)  the  derivatives  of  the  stresses  of 
the  (deterministic)  solution  which  is  computed  by  the  finite  element  method. 
This,  of  course,  needs  a  special  care  and  the  computation  can  be  made  by  the 
postprocessing  technique  (see  [8]). 

The  selection  of  functions  a^,  o^,  does  not  make  usually  any 

problems.  Many  times  we  have  a  traction  free  surface  and  then,  of  course, 
a„  =  =  T_,  =  0  is  the  proper  choice. 

A  y  Ay 

We  have  assumed  that  the  tractions  are  not  stochastic.  We  can  also 
treat  the  combined  case  when  both  the  domain  and  the  traction  are  stochastic. 

We  have  shown  here  only  illustrative  ejcamples  of  relatively  simple 
structure.  The  theory  and  implementation  principles  were  developed  for  the 
general  case.  It  is  possible  to  compute  also  higher  correlation  functions  and 
obtained,  e.g.  the  skewness  of  Che  distribution  of  the  stress  intensity 
factors,  etc. 

In  the  case  of  stochastic  material  coefficients,  we  can  proceed 
similarly  and  get  by  an  iterative  technique  the  desired  data  for  small 
variation  of  the  material  coefficients. 

The  optimal  design  should  in  general  take  into  account  the  stochastic 
character  of  the  input  data. 
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CONCLUSIONS 


Solving  the  problems  of  the  optimal  design  and  the  engineering  problems 
in  general  one  has  to  take  into  account  various  aspects  of  the  mathematical 
model  and  its  numerical  treatment  for  getting  reliable  results.  Detailed  a- 
priori  mathematical  analysis  is  of  utmost  importance  for  the  reliable 
conclusions . 
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The  Laboratory  for  Numerical  analysis  is  an  integral  part  of  the 
Institute  for  Physical  Science  and  Technology  of  the  University  of  Maryland, 
under  the  general  administration  of  the  Director,  Institute  for  Physical 
Science  and  Technology.  It  has  the  following  goals: 

•  To  conduct  research  in  the  mathematical  theory  and  computational 
implementation  of  numerical  analysis  and  related  topics,  with  emphasis 
on  the  numerical  treatment  of  linear  and  nonlinear  differential  equa¬ 
tions  and  problems  in  linear  and  nonlinear  algebra. 

•  To  help  bridge  gaps  between  computational  directions  in  engineering, 
physics,  etc.,  and  those  in  the  mathematical  community. 

•  To  provide  a  limited  consulting  service  in  all  areas  of  numerical 
mathematics  to  the  University  as  a  whole,  and  also  to  government 
agencies  and  industries  in  the  State  of  Maryland  and  the  Washington 
Metropolitan  area. 

•  To  assist  with  the  education  of  numerical  analysts,  especially  at  the 
postdoctoral  level,  in  conjunction  with  the  Interdisciplinary  Applied 
Mathematics  Program  and  the  programs  of  Che  Mathematics  and  Computer 
Science  Departments.  This  includes  active  collaboration  with  govern¬ 
ment  agencies  such  as  the  National  Bureau  of  Standards. 

•  To  be  an  international  center  of  study  and  research  for  foreign 
students  in  numerical  mathematics  who  are  supported  by  foreign  govern¬ 
ments  or  exchange  agencies  (Fulbright,  etc.) 

Further  information  may  be  obtained  from  Professor  1.  Babu^ka,  Chairman, 
Laboratory  for  Numerical  Analysis,  Institute  for  Physical  Science  and 
Technology,  University  of  Maryland,  College  Park,  Maryland  20742. 


